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ABSTRACT 



Aims. To compare simple analytical closure models of turbulent Boussinesq convection for stellar applications with direct three- 
dimensional simulations both in homogeneous and inhomogeneous (bounded) setups. 

Methods. We use simple analytical closure models to compute the fluxes of angular momentum and heat as a function of rotation 
rate measured by the Taylor number. We also investigate cases with varying angles between the angular velocity and gravity vectors, 
corresponding to locating the computational domain at different latitudes ranging from the pole to the equator of the star. We perform 
three-dimensional numerical simulations in the same parameter regimes for comparison. The free parameters appearing in the closure 
models are calibrated by two fit methods using simulation data. Unique determination of the closure parameters is possible only in the 
non-rotating case and when the system is placed at the pole. In the other cases the fit procedures yield somewhat differing results. The 
quality of the closure is tested by substituting the resulting coefficients back into the closure model and comparing with the simulation 
results. 

Results. The simulation data for the Reynolds stress and heat fluxes in the homogeneous case broadly agree with previous compress- 
ible simulations. The closure works fairly well in the homogeneous case with slow rotation but its quality degrades as the rotation rate 
is increased. We find that the closure parameters depend not only on Ta but also on latitude. Furthermore, in the inhomogeneous case 
the closure is unable to reproduce the vertical profiles of the horizontal components of the Reynolds stress. 
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Turbulent convection is responsible for the transport of angular 
momentum and heat in stellar convection zones, in particular in 
that of the Sun. In combination with global rotation, these tur- 
bulent flows lead to the generation of large-scale diff erential ro- 
tation and meridional circulation (e.g. Riidiger 1989), which on 
the other hand, play key roles in sustaining the dynamo of th e 
Sun (e.g. lKrause & Radlei|[!980tlRudiger & Hollerbachll2004l) . 

During the last decades, growing computational resources 
have allowed direct and large-eddy numerical simulations to 
reach a level of sophistication where many aspects of the 
solar differential rotatio n and dynamo can be captured self- 
consi s tentlv (see, e.g. [Miesch et al J 120061: iMiesch & To omre 
2009t iGhizaru et alj|2010t iKapvla et alj|2012l) . However, these 
simulations are still computationally very expensive and cannot 
be employed in performing comprehensive parameter surveys. 
Furthermore, even the currently highest resoluti on simulation s 
are still far from real stars in parameter space (e.g. Kapyla 201 1). 
An alternate way of dealing with the problem is to parametrize 
the small scales by turbulent transport coefficients and solve di- 
rectly only for the large scales. This is often done by approximat- 
ing higher-order correlations by lower order ones in so-called 
closure approaches. In general the pitfall lies in the limited va- 
lidity of the analytical approximations under which the results 
are derived. Hence finetuning of the model parameters is usually 
required. 



Send offprint requests to: J. E. Snellman 
e-mail: jan.snellman@helsinki.fi 



A promising approach to the problem is to validate or cal- 
ibrate the turbulent closure models by comparing their results 
with direct numerical simulations. However, fairly little has been 
done to accomplish this in the astrophysical context. This is es- 
pecially tr ue for closures dea ling with turbulent convection (see, 
however, Garau d et"ai]|2010l) . 

In this paper we build upon previous studies where simple 
analytical closure models were compared with simulations of 
forced or magn etorotationally excited turbulence in fully pe- 
riodic systems (jKiip vla & Br andenbur g 2008; Liliestrom et al.l 
2009: ISnellman et aUl2009t) . Here we extend this work to turbu- 
lent convection in unstratified Rayleig h-Benard setup s , drawi ng 
insights especially from the results of ISnellman et al.l d2012alfbh 
(hereafter S12a and SI 2b) where closure parameters were ex- 
tracted from forced turbulence simulations. Our aim is to com- 
pare three-dimensional direct numerical simulation s (DNS) with 
the clo sure model for convection put forward by iGaraud et al.l 
(1201 Oh (hereafter GOMS10). The bulk of their study is de- 
voted to the derivation and calibration of a closure model for a 
Boussinesq system. Results of DNS and experiments of bounded 
and hence inhomogeneous non-rotating Rayleigh-Benard con- 
vection are referred to for the purpose of determining the free 
parameters of the model. It was found that in the statistically 
stationary state, to a certain extent universal constants can be 
extracted which moreover coincide partly with those from a cor- 
res ponding study of the very different situation of shear flows, 
see IGaraud & Ogilvid (120051) . 

Similarly, an additional free parameter of the closure for ho- 
mogeneous non-rotating Rayleigh-Benard convection was esti- 
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mated on the basis of DNS results where its universality turned 
out to be limited by the destabilization of the fluctuations in 
shallow computational domains. The predictions of the closure 
model for the same setup, but with rotation included, were set 
into relation of previous analytical results. However, no direct 
comparison with corresponding DNS results was performed, in 
particular, there was no independent calibration of the model pa- 
rameters. 

The emergence of coherent structures covering the whole 
vertical extent of the domain was quite generally pointed out 
to be limiting the validity of the essentially local closure. For 
the case of rotating homogeneous Rayleigh-Benard convection a 
further limit was found in the independence of the closure model 
on the rotation rate when gravitation and rotation are perfectly 
aligned. 

Hence, the goal of the present paper is to scrutinize the 
potential of the GOMS10 closure model for rotating convec- 
tion, both for the case of homogeneous and inhomogeneous 
Boussinesq convection. 

2. Models and methods 

2. 1. The Boussinesq system 
2.1.1. Basic equations 

We consider a closure model for Boussinesq convection both in 
a plane layer and infinitely extended space following the pattern 
of iMiller & Garaudl (I2007I) and GOMS10. In general, the time 
evolution of the velocity and temperature fields is governed by 
the Navier-Stokes, continuity, and heat transfer equations 



P^- = - Vp + p(g-2flxU + / visc ), 



= -V-U, 



Dt 
DT 

pcv-^ = V ■ KVT -pV - U + 2vpS 2 , 



(1) 
(2) 
(3) 



where U is the velocity, p is the density, g is the gravitational ac- 
celeration, /vise is the viscous force per mass, p is the pressure, 
T is the temperature, cy is the specific heat at constant volume, 
K is the heat conductivity and D/Dt = d/dt + U ■ V denotes 
the advective derivative. In the following, we assume that g and 
cy have constant values. The viscous force is given by 



U + 2S- Vlnp 



(4) 



with the kinematic viscosity v assumed constant and the trace- 
less rate of strain tensor S, which can be written in component 
form as 

S« = \ (djUi + diUj) - \5 l3 d k U k . (5) 

In the Boussinesq approximation, convection is understood 
as a perturbation to a stationary, purely conductive reference 
state with constant density po which is hence governed by 







Vpo + p g, = xV 2 T + g , 



(6) 



where the thermal diffusivity x = KJ pay is assumed to be con- 
stant and a stationary heat source qo can be included. In this pa- 
per, however, we rely on the simplest case with qo = and a con- 
sequently uniform temperature gradient VXo enforced by appro- 
priate boundary conditions. Denoting the deviations of density, 
pressure and temperature from their reference values caused by 
convection by p' , p 1 and O, respectively, that is, p — po + p', 



p = Po+p' and T = To + O, we obtain from (fJJ and (O for the 
momentum balance 



Po 



DU 

Dt 



= -Vp' + p'g - 2 Po fl xU + p f, 



(7) 



According to the key idea of the Boussinesq approximation, the 
density deviation p 1 from its reference value po is assumed to be 
negligible except in the buoyancy force pg. Density and temper- 
ature perturbations are interconnected by 



L 
Po 



-a(T — To) = — aQ, 



(8) 



where a is the coefficient of thermal volume expansion. Finally, 
the equations of the Boussinesq approximation read 



DU 

Dt 
DO 

~Dt 



= -Vf - a<dg - 2CI x U + v\7 2 U, V • U = 



xv 2 e - u ■ VT . 



(9) 



Here, the reduced pressure $ = p'/po was introduced and the 
viscous force was simplified for the now incompressible flow. 
Contributions to heat transfer from viscous heating and compres- 
sion work are omitted for simplicity. 

2.1 .2. Domain, boundary conditions, control parameters 

For the computational domain we consider a rectangular box 
thought of being cut out at a varying latitude from the convec- 
tion zone of a rotating star. We choose Cartesian coordinates 
(x, y, z) such that their directions locally correspond to those of 
the global spherical coordinates (§, <fi, r) having their axis d = 
aligned with the angular velocity vector O. In (0 the latter has 
then to be written as 



(1 = O (-sintf,0,costf) 5 



(10) 



where t? is the colatitude. We will let the box adopt seven dif- 
ferent positions defined by varying t9 in equidistant steps of 15 
degrees from 0° (pole) to 90° (equator). For the dimensions of 
the box, L x , L y , L z , we assume L x = L y while L z may differ 
from them. 

For all quantities, periodic boundary conditions in the hori- 
zontal (x, y) directions are employed throughout the paper. For 
the conditions at the z boundaries we adopt a combination of 
non-penetrating stress-free flow/zero temperature perturbation, 
i.e. 



d z U x =8 z U v = u z = e = o at 



0,L Z (11) 



or periodic conditions for all quantities. While in the former case 
the constant reference temperature gradient VXb is given by the 
reference temperature difference across the layer, it has to be 
considered an independent parameter of the equations ^ in the 
latter. The non-periodic case shows the complication of the con- 
vective turbulence becoming inhomogeneous near the bound- 
aries. In contrast, the periodic case does not feature any causes 
for inhomogeneity and is therefore also labelled as homogeneous 
Rayleigh-Benard convection. In both cases the turbulence natu- 
rally has to exhibit anisotropy already because of the preference 
of the gravity, i.e., the z direction. 

The system (0 is governed by the following three dimen- 
sionless parameters: The magnitude of the temperature gradient 
(and eventually the vigour of the convection) is quantified by the 
Rayleigh number 

Ra=^ 3AT °, (12) 
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where g = g z , ATo is the reference temperature difference be- 
tween the top and the bottom of the domain and d = L z its ver- 
tical extent. For the case of periodicity in z, the definition ( fT2] > 
has to be modified properly employing the constant background 
temperature gradient V z Tq = Go, that is 



Ra 



agd i Go 
vx 



(13) 



The ratio of viscosity and thermal diffusivity is given by the 
Prandtl number 

(14) 



Pr = - 

X 



and finally the rotation rate is measured by the Taylor number 



(15) 



Another way to express the strengths of rotation and viscous ef- 
fects, but in the form of diagnostics rather than of control param- 
eters, is provided by the Coriolis and Reynolds numbers, respec- 
tively, 



Co = 



Re = 



vkf 



where k{ = 2ir/d is the wave number corresponding to the ver- 
tical extent d and U Tras = (J„ U 2 dV jV) 1 ! 2 is the root mean 
square velocity with the volume of the domain V. 

2.1 .3. Closure model 

In deriving the closure model we first have to specify the 
averaging procedure by which the mean quantities are de- 
fined. Given the periodicity of our models in the horizontal 
directions, averaging over these directions is appropriate and 
we will apply it throughout this paper. Hence, the mean of 
a quantity /, indicated by an overbar, is given by f(z) = 
II II f{ x TVi z )dxdy/L x Ly. This procedure satisfies all the 
Reynolds averaging rules. Denoting turbulent quantities with 
lowercase letters, we have U = U + u, = 9 + 6* etc. and 
the equations for the mean quantities read 



V • U 

dU 

dt 



= 0, 



(17) 



-U ■ VU - aOg - V* - 2fl x U + vV 2 U 
-VIZ, (18) 



de 
dt 



— = -[/• V(T + 6) + X V 2 9 - V • T, 



where 7Z is the Reynolds stress tensor, IZij = Uiiij, and T is the 
turbulent heat flux, Ti = 9ui. These equations are identical to 
those used in GOMS10, with the exception of the extra Coriolis 
term. As a consequence of the chosen average, horizontal_deriva- 
tives vanish and the continuity equation reduces to d z U z = 0, 
from which it follows that U z = const. For the plane layer with 
impenetrable boundaries this means U z = 0. The equations for 
the remaining components of the mean velocity read, with grav- 
ity having only a z component, 

t x = -d z K xz + 2fl z U y + vd zz U Xl (20) 

t y = -d z n yz + 2(n x u z - n z u x ) + vd zz u y . 



Note that we do not need to solve for the_mean pressure ^ as it 
would only affect U z . The equation for 9 reduces to 



9 = X d zz Q - d z F z 



(21) 



Evolution equations for the Reynolds stress and turbulent heat 
flux can be derived from the equations for the turbulent quanti- 
ties u and 9, see Appendix lAl In doing so one comes inevitably 
across higher-order correlations of m and 9. The essential step of 
the closure procedure as proposed in GOMS10 consists then in 
replacing these correlations by aggregates of second-order corre- 
lations, more specifically, of the quantities 7Z and T themselves. 
In an analogous way, some second-order correlations which can- 
not directly be expressed by the components of 1Z and J~, or Q 
are modelled. The details can again be found in Appendix lAl In 
their final form the closed set of equations reads 



Hi 



U k d k Kij + UikdkU j + TZjkdkUi + affigj + F.jg. t ) 
vdkkTtij + 2Qi(euk'R-jk + £jikHik) (22) 



C l + C 2 
L 



n l/2 



3L 



(16) ti + UjdjFi + TidjUi + Kijdj (9 + T ) + aQ 9i 



(23) 



a 



L 2 



where 1Z = 1Z X 



K zz is the trace of 71, Q = 9 2 



the temperature variance, C\ t 2,6,7,u,vx,x are model parame- 
ters and L is a characteristic length scale, here taken to be the 
distance to the closest boundary, L = min(z, d — z). Hence, 
the closure terms are explicitly position dependent. In order to 
obtain a closed system of equations, (f20l> — (f23T> have to be com- 
pleted by an additional one for Q which can be derived using 
Eq. (0. With a closure applied in the same way as for TZ and 
T, it can be written as 



Q + U i d i Q + 2F i d i {Q + T ) - xduQ 



L 



(24) 



As for the original problem, we assume stress-free impene- 
trable vertical boundaries at z = 0, d for the mean velocity and 
given that corresponding conditions hold for the turbulent veloc- 
ity u, we may require 



(19) d z 1Z xx = d z 1Z x 



dzTZ-u 



1Z XZ = 1Z yz = 1Z ZZ = 0, (25) 



although the first three are only necessary, but not sufficient. 
Given that the temperature T_is fixed at the boundaries, the per- 
turbation 9 and with it both 9 and 9 must vanish there. Thus the 
boundary conditions for mean temperature perturbation, turbu- 
lent heat flux and temperature variance are 



Q = J r x=J 7 y=J 7 z = Q = 0- 



(26) 



For homogeneous turbulence which is obtained in the fully 
periodic setup, the closure model can be further simplified by 
defining mean quantities as volume averages over all space or, 
what is the same, over the computational domain. Then the 
second-order correlations TZij, Ti, Q are completely decoupled 
from the mean fields U and 9 and for the latter only trivial 
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solutions exist, see (f20b and (f2TT) . The remaining set of ordi- 
nary differential equations with respect to time is provided in 
Appendix |U see also GOMS10, Sect. 4.4. Note that now there 
is no longer an argument to identify the length scale L in the 
closure ansatzes with the distance to the nearest boundary, as the 
periodic boundaries do not limit structure formation. We have 
chosen to identify L with the vertical extent of the computational 
box (the fundamental period of the considered fields in z ), but it 
could also be considered a free parameter. Anyway, by compari- 
son with DNS or experiments only the combinations Cifl&^/L 
and C Vy „ xa /L 2 can be determined. 

2.2. DNS setups 

For the DNS two different setups are used, both in a local 
Cartesian volume of size L x fully periodic setup 

with a uniform background temperature gradient as described in 
Sec. l2.1l for the homogeneous case, and an inhomogeneous setup 
where either stress-free 

d z U x = d z U y = U z = 0, (27) 
or no-slip 

U x = U y = U Z = 0, (28) 

boundary conditions are applied for the velocity. The tempera- 
ture is kept fixed at the vertical boundaries. In the inhomoge- 
neous setup the aspect ratio of the box is five, i.e. L x = L y = 
5L Z , whereas it is one in the homogeneous setup. 

The numerical simulations were performed with the PENCIL 
which uses sixth-order accurate finite differences in 
space, and a third-order accurate time-stepping scheme, see 
iBrandenburg & Doblerld2002l) ; lBrandenburg|(l2003l) . In all cases 
the time integration was advanced until a statistically station- 
ary state was reached. In addition to the horizontal or volume 
averages, a time average over this state is taken because in par- 
ticular in the homogeneous case the spatial averages still show 
strong fluctuations. Errors are estimated by dividing the time se- 
ries into three equally long parts and computing mean values 
for each part individually. The largest departure from the mean 
value computed for the whole time series is taken to represent 
the error. 

3. Results 

GOMS10 stated strong rotation to be an effect that could pos- 
sibly reduce or completely rule out the applicability of the clo- 
sure. The onset of convection and its saturated stage as a func- 
tion of rotation were studied by Mille r & Garaudl d2007l) using 
essentially the same closure model as that in GOMS10, but they 
did not directly compare mean quantities like TZ from the model 
with the corresponding ones from DNS. Hence it is the main pur- 
pose of this paper to study the behavior of the GOMS10 model 
for a range of rotation rates (or Taylor numbers) and for differ- 
ent inclinations of the rotation axis by performing such compar- 
isons. 

3.1. DNS Runs 

The homogeneous runs are summarised in Table Q] for more de- 
tails see Table IB. II in the Online Material. In each of the sets 

1 http://code.google.eom/p/pencil-code/ 



Table 1: Summary of the homogeneous DNS runs at different 
colatitudes Ra = 3 • 10 5 and Pr = 1 throughout. Ta = 
4 • 10 4 . . . 1.3 • 10 7 in sets A through G . 



Set Re Co U rms /(dagAT ) 1/2 



z 






91 







0.57 


A 


0° 


91 


- 127 


0.06- 


0.87 


0.57-0.80 


B 


15° 


81 


-91 


0.06 - 


1.09 


0.51 -0.57 


C 


30° 


70 


-89 


0.06 - 


1.07 


0.44 - 0.56 


D 


45° 


62 


-91 


0.06 - 


1.39 


0.39-0.57 


E 


60° 


57 


- 96 


0.06 - 


0.95 


0.36-0.61 


F 


75° 


54 


- 89 


0.06 - 


1.02 


0.34 - 0.54 


G 


90° 


48 


- 85 


0.06 - 


1.22 


0.30 - 0.54 



A-G the latitude was kept fixed, but the rotation rate was var- 
ied. Z denotes the non-rotating run. The ranges for the Reynolds 
and Coriolis numbers as well as the rms velocity U rms probed 
by DNS are listed. 

In the homogeneous DNS runs we had to deal with a nu- 
merical stability problem in the transition from the kinematic, 
exponentially growing, stage to the stationary stage if the Taylor 
number was too small. To circumvent this we started the runs 
with high (« 10 7 ) values of Ta which allowed a statistically sta- 
tionary state to be established. Then we gradually lowered Ta 
until the desired parameter range had been reached. Eventually 
we were able to successfully perform even a non-rotating run us- 
ing this method. A snapshot of the vertical velocity in this Run Z 
is pictured in the left panel of Fig. Q] In the time series of the 
statistically stationary state, we observe large fluctuations and 
intermittent exponential growth, most likely a s a manifestation 
of the so-called "elevator modes", described in lCalzavarini et all 
(2006), which are weakly damped for a box aspect ratio of unity 

As a more realistic example, we performed simulations with 
an inhomogeneous setup with the similar values of Ra and Pr as 
in the homogeneous runs. The distribution of the velocity com- 
ponent U z of a non-rotating DNS run with stress-free boundary 
conditions is shown in Fig.[T| right panel. One can immediately 
see that there are flow patterns encompassing the entire vertical 
span of the simulation box. As the closure ansatzes in (I22l-(l24li 
are local in space and time, the presence of such large-scale co- 
herent structures can adversely affect their applicability, as al- 
ready noted by GOMS 10. 

3.2. Homogeneous case 

3.2.1 . Calibration of the closure model 

We assume that the Rayleigh numbers are high enough so that 
we can set C v = C vx = C x = throughout. This is prob- 
ably justified in the case of the homogeneous setup where no 
boundary layers are present although a parameter scan in Ra 
would be required to ascertain this for sure. In the inhomoge- 
neous case this is less clear. However, apart from boundary lay- 
ers there are other more serious issues that affect the results and 
that will be described in Sect. I3.3I Stationary solutions of the 
closure model, Eqs. ( IB. II ) in Appendix iBl result for given {Cj} 
from the corresponding nonlinear algebraic system of equations. 
This opens up systematic ways of calibrating the parameters. 
Two such methods are described here. 
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Run Z ■ ni 41 



Fig. 1: Velocity component U z in units of (dagATo) 1 / 2 from non-rotating DNS runs. Left: homogeneous case Run Z, Ra = 3 • 10 5 ; 
right: inhomogeneous case Run Zsf, with stress-free vertical boundaries, Ra = 5 • 10 5 . 



Least squares fit By inserting the parameters (ag, f^o, Go) 
used in the DNS runs together with their statistically stationary 
results for TZ, T and Q into the stationary version of the clo- 
sure equations dB.lt . a generally overdetermined system of lin- 
ear equations for the Ci.2,6,7 ( ten equations vs. four variables) 
of the form 

Nc = L (29) 

is obtained where c = (C±, C2, C$, C7). The matrix TV can 
be derived from the closure terms and the vector L contains 
all remaining terms, such as the Coriolis and buoyancy terms, 
see Eqs. ( IB.lt . At the pole, however, the system d29l is no 
longer overdetermined and c can be determined unambigously 
as demonstrated in Appendix IB"1 At all other latitudes this sys- 
tem can be solved only approximately using the standard lin- 
ear least squares method, that is, solving the regular system 
N T Nc = N T L, where the superscript "T" denotes transpo- 
sition. The results from this procedure are shown in Fig. [2] as 
functions of Taylor number and colatitude. As an indicator of 
the quality of the solution, we also calculated the residual norm 

Rest = \\Nci s -L\\, (30) 

where the subscript "Is" refers to the least-squares solution and 
1 1 • 1 1 denotes the Euclidian norm, see Fig. [2] Another way to 
quantify the quality of the obtained closure model consists in 
calculating the difference in the stationary solutions for X = 
(7l X x, . . . , TZ zz ,F x ,F y , T z , Q) from the closure model dB.lt 
with c = c\ s and from the corresponding temporally averaged 
results of the DNS run, that is, calculating the residual 

Resx = 1 1 -^closure — -XdnsII- (31) 

This quantity is also shown in Fig. [2] together with the left hand 
side of the realizability condition 

2C 6 - C 7 - Ci - C 2 > (32) 

given in GOMS10 and the largest eigenvalue of the linearized 
system, see below. The solutions X c i osme and Xdns are directly 
compared in Fig. [3] When Eq. d29l ) is overdetermined, there is in 
general not any set of coefficients {Ci} with which the closure 
reproduces all the modeled quantities perfectly. This is indicated 
by the residual (fJTJ. 

From Fig. [2] one can see that the derived model coefficients 
change with colatitude and Taylor number, with following pat- 
terns: at the pole they fall with growing Ta, but not so for any 
other colatitude. Instead, they first grow with Ta and plateau or 



Table 2: Ratios of the coefficients Ci obtained from the non- 
rotating Run Z compared to GOMS10 and S12a,b. 



Ci/Cj 


this study 


S12a 


S12b 


GOMS10 


Ci/Ca 


0.82 


w 0.4 




0.66 


Ci/Cg 


0.42 






0.29 


Ci/C, 


0.58 






0.29 


C 2 /C 6 


0.51 






0.43 


C2/C7 


0.72 






0.43 


C 6 /C 7 


1.40 




as 2.1 


1.00 



fall for Ta > 10 6 . Both growth and fall get steeper with growing 
colatitude, and all the curves converge as Ta approaches zero. 

Using the (exact) results for the non-rotating run we com- 
puted the different ratios of the coefficients {Ci} and compared 
them to the corresponding ratios from GOMS 10, S 12a, and S 12b 
in Table 12 These ratios are important because any difference in 
the non-rotating case in the results for C, could be due to a badly 
chosen lengthscale L, which is canceled by the ratios. In any 
case, we see that our values are at odds with those of GOMS 10. 
As for the residuals, they unsurprisingly vanish (up to round- 
off errors) at the pole, otherwise rise with Taylor number while 
the colatitude has only a small effect on the residual (l30l l. but a 
stronger one on (l3Tt . At the equator the residual Resx settles at 
unity for Ta > 5 • 10 7 . This is because the analytical results from 
the closure become very small for all the modeled quantities in- 
dicating that the obtained {Ci} do not allow for another than the 
trivial solution of dB.lt . 

The DNS results Xdns for a slow and a rapid rotation case 
are compared with the corresponding closure results obtained 
with the derived model coefficients in Fig. [3] At slow rotation the 
DNS and the closure results are visibly closer to each other than 
at fast rotation. Further, the fit for small quantities like 7Z X y,xz,yz 
or T x>v is significantly worse than for the large ones, like 7Z ZZ , 
T z or Q. This is mainly because residuals of large quantities 
contribute stronger to the norm of the residual vector which is 
minimized by the least squares method. 

Optimization approach In order to determine the {C;} from 
the results of the DNS in an optimum way, we can also formu- 
late the following optimization problem: minimize the objective 
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function 

(-X"closurc — -XdNs) (33) 

obeying the constraints Ci,2,6,7 > and d32b . The prob- 
lem was solve d by th e Generalized Reduced Gradient Method 
(lLasdon et alJ (1 1978h ~) as implemented by the IDL rou- 
tine CONSTRAINED.MIN with the nonlinear system from 
Eqs. ( IB.U being solved by Newton iteration (IDL routine 



NEWTON). The optimum results are shown in Fig. [4] where the 
upper four panels refer to the {Ci } while the lower left one gives 
the normalized objective fucntion 

(-^closure — -^DNs) 2 /-^DNS (34) 

and the lower right one the value of the quantity 2Cq — C7—C1 — 
C2 from the constraint d32b . Along with the dependence on Ta, 
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Fig. 3: Results of the closure model with coefficients {C, } from the least squares fit compared with the corresponding DNS results. 
Note that the {Ci} depend on rotation rate fig and colatitude see Fig. [2] Dotted lines/squares: closure/DNS results for slow 
rotation (Runs A2-G2 with Ta = 1.6 • 10 5 ), solid lines/diamonds: for faster rotation (Runs A7-G7 with Ta = 1.96 • 10 6 ). 



there is again in general a separate one on -d. At the pole (z9 = 0) 
the objective function assumes exceptionally low values. This is 
a consequence of the already mentioned degeneration which al- 
lows to determine the {Ci} uniquely from the Xdns- Hence the 
objective function actually vanishes and the observed values are 
completely due to the iterative nature of the solution and round- 
off errors. Apart from the pole, the quality of the optimum is 
in general decreasing with growing Ta, yet having only a weak 
dependence on colatitude. For 10 6 < Ta < 1.6 • 10 7 and Ta 
dependent i? intervals between 45° and 75°, the constraint (l32l 
becomes "active" in the sense that the optimum lies then on the 
margin of the admissible domain, that is, 2Cq — C\ — C^ — CV = 
(red curve sections in Fig.|4]i. 

Given that, apart from the constraint (l32l , also the stability 
properties of the closure model should not differ from that of 
the DNS, we enhanced the optomization problem by the con- 
straint that for the optimum fit the stationary solution corre- 
sponding to it should be stable. For that, we linearized the sys- 
tem (IB.U about the state X c i osurc , obtaining a system of the 
form dt (SX) = A ■ SX for the perturbations dX, and required 
that the maximum of the real parts of the eigenvalues of A is 
negative. To avoid influences of numerical noise we set their up- 
per bound to a small negative value instead of zero. The matrix 
eigenvalue problem was solved by means of the IDL routines 



LA_ELMHES and LAJHQR. It turned out that the additional 
constraint never becomes active, that is, that stability is already 
granted if only d32b is obeyed. 

In Fig. [5] the DNS and closure model results as functions 
of the colatitude are given for two different Taylor numbers. 
Obviously, the dominant variables 1Z ZZ , J- z and Q are always 
acceptably fitted. For slow rotation, also the other quantities ex- 
cept lZ yz show good fits. By applying appropriate weights in the 
objective function d33l the fit quality could be distributed "more 
equally", but such a procedure suffers from high arbitrariness. 

Both versions of the optimization approach produce up to 
roundoff errors identical results as long as neither of the two 
constraints is active. The dependence of the objective function 
on latitude is in general weak, and beyond Ta ps 10 7 its depen- 
dence on Ta is weak too. Again, the ratios C1/C2 and C^jOi 
for the non-rotating case are at odds with GOMS10, see Table[3] 
When comparing the results of the least-squares and the opti- 
mization approaches we find major quantitative differences in 
Ci,2,7> m addition differing monotony in C1.7 while the residu- 
als are clearly smaller for the latter approach, namely ps 30% vs. 
up tp 100 % for the former. 

Dependence on Ta For both approaches, the obtained {C. t } 
show a strong dependence on Ta. At the pole, however, there is 
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Table 3: Ratios of the coefficients Ci obtained for rotating runs 
by the optimization approach. Ta = 4 • 10 4 . . . 10 8 



■& 


C1/C2 


Ci/Ce 


C1/C7 





0.81 


- 1.73 


0.42 


-0.52 


0.58 


- 1.45 


15 


0.55 


-0.88 


0.29 


- 0.42 


0.55 


-0.88 


30 


0.52 


-0.81 


0.21 


- 0.42 


0.52 


-0.81 


45 


0.54 


-0.83 


0.13 


- 0.42 


0.51 


-0.83 


60 


0.51 


-0.89 


0.06 


- 0.42 


0.51 


-0.89 


75 


0.49 


-0.89 


0.01 


- 0.42 


0.49 


-0.89 


90 


0.43 


- 1.11 


0.01 


- 0.42 


0.43 


- 1.11 



ansatz would also introduce a dependence on the rotation rate f2rj 
(or Ta). Possible additional closure terms are, e.g. 

r^rij, l^n^liVtj -\- T^ji^li^li for T^ij j 
Eijk^ljJ-k for Ti , 

where all corresponding coefficients may depend on \g ■ fi| 
and/or |fi|. 

3.2.2. Reynolds stress and heat fluxes in comparison to 
compressible simulations 



The off-diagonal Reynolds stresses and turbulent heat flux are 
no such explicit dependence in the solution ( IB.1U . We interprete important in generating the differentia l rotation of stellar con- 
this conflict as evidence that the isotropic ansatz for 1Zij in d22l > vective envelopes (e.g. Rudiger 1989). These quantities have 
is quite generally inadequate in view of the anisotropy induced been computed from numerou s simulations of com p ressible con- 
by the direction of the rotation vector. An appropriate anisotropic vection in Cartesian (e.g. Pulkk inen et alJ fl993t IChanl 1200 lb 
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Fig. 5: Results of the closure model with coefficients {C,} from the optimization approach with the corresponding DNS results. 
Dashed lines/squares: slow rotation (Ta = 4 • 10 4 ), solid lines/diamonds: faster rotation (Ta = 1.3 • 10 7 ). 



Kaovl a et al.ll2004t iRtidiger et alJfe OOSb) an d sphe rical geome- 
tries (e.g. lRieutord et aO 19941: iKapvla et alj|201 lb . It is impor- 
tant to compare the results of our homogeneous Boussinesq runs 
to those in the literature in order to draw conclusions on the ro- 
bustness of certain features such as the latitude and rotation rate 
dependence. 

We find that lZ xy , corresponding to latitudinal flux of an- 
gular momentum is always positive, i.e. directed towards the 
equa tor in accordance with previous DNS and analytical the - 
ory (Kichatinov & Riidiger 1993; Kitchatinov & Rudiger 2005). 
There is a tendency for the maximum of lZ xy to move toward the 
equator as the rotation rate is increased in accordance with com- 



pressible simulations of lChanl (l200lh and iKapvUTeTa l. (2004). 
Furthermore, the vertical flux corresponding to lZ yz is always 
negativ e. No sign reversal, observed at high Co in compressible 
runs of iKanvlaet all (12004 . is seen even for the highest Taylor 
numbers. The third off-diagonal component 1Z XZ is mostly neg- 
ative, although positive values occur at mid-latitudes for rapid 
rotation. Earlier re sults suggest that posi tive values occur at high 
latitudes only (e.g. Pulkkine n et al.ll 19*931) . 

The latitudinal heat flux T x is always directed towards the 
pole and the azimuthal heat flux T v is negative, i.e. in the ret- 
rograde direction. These features are broadly in accordance with 
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Fig. 6: Vertical profiles of the stationary solution for the inhomogeneous non-rotating case with stress free (solid lines) and no slip 
(dash-dotted) boundary conditions. Dashed curves: corresponding results from a representative ID closure model with stress-free 
boundaries. 



Carte sian (e.g. Rudiger et al.ll2005ah and spherical simulations 
(e.g.|KapylaetaJj2011|). 



3.3. Inhomogeneous case 

Our numerical experiments with the inhomogenous DNS setup 
included rotating runs at different latitudes as well as non- 
rotating runs, though only the latter will be considered in this pa- 
per. The main reason for excluding the former is that the closure 
and DNS results could not be brought to a satisfactory agreement 
already in the non-rotating case. The nature of this problem is in- 
dependent of rotation and will be elaborated in the next Section. 

For a good calibration of the closure model we have to re- 
quire that in the statistically stationary state the vertical (z) pro- 
files of the mean quantities from the DNS run under considera- 
tion and those from the stationary solution of the corresponding 
closure model should be in decent agreement. Thus, Pr and Ra 
(see Eqs. (fTZl) and (O) were set equal in both setups, an initial 
choice for the {C\} (with C VtK>VK neglected) was made and the 
system d20l)-(l24l). discretized with second order difference for- 
mulae, was Euler integrated in time until a stationary state was 
reached. Then the z profiles were compared to those of the corre- 
sponding DNS run. As a sanity test we checked in each timestep 
whether the diagonal components of the Reynolds stress tensor 
remained positive, and terminated the run if not. 

We show the results for all non-vanishing quantities, i.e. 
IZyy, 1Z ZZ , T z , Q and 9 for two non-rotating runs (labeled 
Zsf and Zns) with Ra = 5 • 10 5 and Pr = 1 in Fig.0 In Run Zsf, 
stress-free boundary conditions ( |27T ) are applied for the velocity 
whereas in Run Zns, no-slip boundaries d28l ) are used. We per- 



formed a representative ID closure model run where the ratios 
of the Ci parameters were taken from the non-rotating homoge- 
neous run (see Table |2]i while we scaled C\ so that fair agree- 
ment for F z between the DNS and closure is obtained. Stress- 
free boundaries are applied in the closure model. The resulting 
DNS vs. closure model profiles for the non-zero quantities are 
also shown in Fig. [6j w_herefrom we see that the general forms 
and magnitudes of JF Z , 8 and 1Z ZZ are in qualitative agreement, 
while the shape of the profile of Q is qualitatively correct but 
the magnitude from the closure model is much smaller than that 
from DNS. The most severe discrepancy is observed for the di- 
agonal components 7Z XX and 7Z yy , for which neither the profiles 
nor magnitudes are correctly reproduced: The DNS results form 
cuplike (U) profiles, while the profiles from the closure model 
resemble more caps (n). Although we have set C„, C x and C vx 
to zero, the basic picture stays the same even when these pa- 
rameters are given non-zero values, likewise when changing the 
other parameters. 

For stress-free boundary conditions in the DNS models the 
coexistence of cuplike results for 7Z. XX and lZ yy with caplike re- 
sults for 1Z ZZ is physically plausible: 1Z ZZ vanishes at the bound- 
aries because the flow cannot permeate. When descending flows 
hit the bottom or ascending flows hit the top, they are deflected 
and redistribute their kinetic energy from the vertical compo- 
nent to the horizontal ones which need not vanish at the bound- 
ary. Hence one would expect 1Z XX and lZ yy to rise when 1Z ZZ 
tends to zero near the boundaries. The situation can potentially 
be improved by using no-slip boundary conditions for the veloc- 
ity, forcing the horizontal velocities to obtain zero values at the 
vertical boundaries; results from such a DNS run are shown in 



10 



Snellman et al.: Testing turbulent closure models with convection simulations 



Fig.|6]with dash-dotted lines. This has the desired effect near the 
vertical boundaries, but 7Z XX and lZ yy retain their cuplike pro- 
files away from the boundaries. 

The mathematical reason why 7Z XX and lZ yy from the 
closure model decrease near the boundaries is related to the 
nature of its basic ingredients: the relaxation terms of the 
form C\R}l 2 L~ 1 'R.ij and the isotropization terms of the form 
C 2 n 1 / 2 L- 1 (n i , J - \5ij1l). Without any other influences, the 
relaxation terms tend to cause the stresses to vanish, while 
the isotropization term forces the diagonal components of the 
Reynolds stress to have the same profile: Let us examine 1Z XX 
in the stationary case ignoring the viscous terms for simplicity. 
Since U vanishes in this case, we get from Eq. (I A. 1 Oi l 

^=krc^ (35) 

independent of the prescription of L(z). So 1Z XX simply has the 
same profile as 1Z and for lZ yy the situation is exactly the same, 
making them equal and both proportional to 1Z ZZ . Hence, all 
three stresses have to vanish at the boundary if only 1Z ZZ has to 
do so, that is, at any impenetrable boundary. When including the 
closure term with C v , 1Z XX and lZ yy are no longer proportional 
to 1Z ZZ , but still vanish with it. A mechanism to transfer kinetic 
energy between diagonal components of the Reynolds stress ten- 
sor near the boundaries is missing in the GOMS10 model what 
becomes obvious if they are stress-free. 

From Eq. (38) of GOMS10 it can even generally be con- 
cluded that, independent of the boundary condition for U, any 
regular solution of the closure model must exhibit vanishing 1Z at 
the boundaries as long as the length L is defined as the distance 
to the closest boundary. 

4. Conclusions 

The main conclusion to be drawn from our experiences with the 
homogeneous Boussinesq model is that the closure parameters 
{Ci} vary as functions of Taylor number and latitude what can- 
not be captured by a dependence on a single quantity, say, | SI ■ g \ . 
Further, the validity of the closure degrades as the rotation rate 
and colatitude are increased, as indicated by the growing resid- 
uals of the parameter fits (see Figs. [2] and |4}. Both facts sug- 
gest that the GOMS10 closure is essentially incomplete in the 
rotating case. In particular, the purely isotropic (even isotropiz- 
ing) character of the closure ansatz should be revised, given the 
clear anisotropies which are induced by the direction of rotation. 
However, even in the non-rotating case we were not able to re- 
produce the ratios of the coefficients {C\} provided in GOMS 10, 
see Table [2] This might be tracked down to the fact that in 
their study the coefficients were not independently determined 
from DNS, but adopted without change from the inhomogeneous 
model while only adjusting the length scale L. 

We observed that positivity of the parameters and their ad- 
herence to the realizability condition alone always guarantees 
stability of the stationary solutions of the closure model in ac- 
cordance to the stability of the underlying statistically stationary 
DNS solution. 

In general we have not studied Rayleigh numbers higher than 
Ra ~ 5 • 10 5 due to computational constraints that arise as a con- 
sequence of the higher resolution required. This might be rele- 
vant for the comparison to GOMS 10 having employed up to two 
orders of magnitude higher values. Another aspect of potential 
importance for such comparison is our computational box as- 
pect ratio of unity (they have less than unity) which nevertheless 
allows for statistically stationary solutions in DNS. 



For the inhomogeneous setups, too, these findings keep their 
validity: we could not reproduce the GOMS 10 results quantita- 
tively, but our results are again limited by the moderate Rayleigh 
numbers. The model parameters were found to be not univer- 
sal, but dependent on rotation rate and latitude. However, we 
abstained from presenting the latter findings in the paper as al- 
ready in the non-rotating case we detected a basic difficulty: it 
was virtually impossible to choose model parameters such, that 
the profiles of 7Z XX and 7Z yy from the DNS could be reproduced 
qualitatively by the closure model although this was success- 
fully demonstrated in GOMS 10. As a reason for this discrep- 
ancy we could identify the stress-free boundary conditions of 
our model as opposed to the no-slip ones in GOMS 10: by virtue 
of its isotropization terms, the closure ansatz is unable to model 
the redistribution of kinetic energy from the normal to the lateral 
flows close to any impenentrable boundary. This flaw is veiled 
for no-slip conditions because all flow components then tend to 
zero near the boundary, but not so for stress-free conditions. As 
a consequence, we obtained convex instead of concave z pro- 
files for the mentioned Reynolds stress components being just 
those which do not vanish at stress-free boundaries. As a way 
out, terms should be added to the closure which explicitly model 
the named energy redistribution, that is, which are of explicit 
"anisotropizing" nature and are in effect only close to the bound- 
aries. Most likely the closure model could profit from such terms 
for any impenetrable boundary condition. 
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Appendix A: Closure model equations for the 
inhomogeneous Boussinesq system 

The equations for the fluctuating quantities u and 9 read 



du 
~dt 

dt 



- U ■ Vu - u ■ VU - a9g - Vtjj 

- 2f2 x u + vV 2 u - V ■ (u <E> u - 71), (A.l) 
-U-VO-u- V(6 + T ) + xV 2 - V • (9u - T), 



with ® denoting the dyadic product, or, in component form 



~dt 

00 
di 



UjdjUi — UjdjUi — aOgi — diip 

2£ ijk tt j u k + vdjjin - dj(uiUj - Kij), (A.2) 



= - u^t 



ujd^e + To) + v) u (> - <)■(>,, , - T j). 



Apart from the Coriolis term, these equations are identical to 
those given in GOMS 10. From them we can derive the following 
equations for the Reynolds stresses IZij, the turbulent heat flux 
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T% and the temperature variance Q (with U z = 0) 



Tty + U kdkTZij + HikdkU j + U jk d k U l + affigj + T^gC) 
— vdkkT^ij + 2Qi(eukTijk + £jikR>ik) = (A. 3) 



- Uidjip + Ujdiip - Uid k {ujU k ) + Ujd k (u t u k ) - 2vd k u i d k u jl 
- \{v + X)djjFi + 2e l]k Vl J F k = 



(A.4) 



— Qdi\j) — 8dj{uiUj) + Uidj(6uj) 

+ \( v ~ x)dk{0d k u l - Uid k 6) -(v + x)d k 8d k u il 

Q + U l d l Q + 2Fidi (9 + T ) - X da Q = (A.5) 



2ed i {6u i )-2 X {d i 6y. 



Here, the right hand sides contain third order correlations of tur- 
bulent quantities (including the correlations with the pressure iff) 
and terms originating from the Laplacians which cannot be ex- 
pressed by the considered second order correlations. In the clo- 
sure model of GOMS10 all these are replaced in the following 
way: 



Uidjil) + Ujdiip + Uidk(ujU k ) + Ujdk{uiUk) + 2vd k u l d k Uj 



= A lz n ij - < £n 3 / 2 s ij , 

0±J 



(A.6) 



9dii/j + 9dj(uiUj) + Uidj(9uj) 

-\{v - x)d k {0d k ui - Uidf.9) + {v + x)d k 0d k u l 



n 1 ' 2 F l + Uu + x)^P = ^P, 



L 



L 2 



(A.7) 



2ed i (eu i ) + 2 X {d i ey -> ^-h 1 / 2 q+ x ^q = a q qca.8) 



with 



A 



_ (Ci + C 2 ) 1/2 c v 
n- z ft + VJ? , 



A r = ^TjVa + i (l/ + x) ^ ; Ag = + x 11. 



(A.9) 

L 2 



T^xz 

V — 

V — 

n zz = 



-A n K xx + §^ 3/2 , (A. 10) 

^-xz^zU y T^-yz^zU x ~t~ ^^x^-xz 

+2n z (TZ yy - K xx ) + vd zz K xy - A R TZ xy , (A.ll) 

^-zz^zU x OtJ~ xQz ^^x^-xy H~ ^^z^yz 

+vd zz 1l xz - A n K xz , (A. 12) 

'ffl^yz^zUy ~\~ 4:£l x 7^yz A£l z 7^, X y H~ ^^zz'^yy 

-A n n yy + §^ 3 / 2 , (A. 13) 



TlzzdJJy - aT y g z + 2Q X (1Z ZZ - 72^) - 2fi 2 7\L a;z 
vd zz TZ yz ~ A n Tl yzl (A. 14) 

2aF z g z - iCl x TZ yz + vd zz 1l zz 



(A. 15) 



T x 



T z = 



-K xz d z (0 + T ) - T z d z U x + \{v + X )d zz T x 
-2Sl z T y -A T T x , (A. 16) 

-K yz d z (e + T ) - T z d z U y + \{v + x)dzzf y 
-2£l x F z - 2Vt z F x - Aj^Fy, (A.17) 

-Tl zz d z (Q + T ) - aQg z + \{v + X )d zz T z 
-2n x F y -A T F z , (A.18) 



Q = -2F z d z (Q + T ) + xd zz Q-A Q Q., 



(A. 19) 



Given that the temperature profile To is linear in z and_ steady, 
the last equations andjhe evolution equation dT9b for 9 can be 
simplified by setting Q + Tq — > 9. Then the boundary condi- 
itons for 9 of have of course to be changed to the original ones 
for T. Note that the stationary version of the autonomous system 
(IA.10t - dA.19t with (l20"li-(l2TTi does have non-trivial solutions as 
demonstrated in GOMS10. Due to the nonlinearity of the sys- 
tem they exist not only for specific combinations of its parame- 
ters like in linear eigenvalue problems, but (at least within wide 
margins) for any specification of them. 



Thus, the closure consists of relaxation terms, such 
as C 1 L- 1 Tl 1/2 'Ri j , isotropization terms C 2 L- 1 K 1/2 (R. tj - 
^IZSij) and terms like vC u L~ 2 Hij corresponding with turbu- 
lent diffusion. For the length scale L we adopt the distance to 
the closest vertical boundary. Applying these simplifications we 
arrive at the equations (l22l424l i. In component form they read 



Appendix B: Closure model equations for the 
homogeneous Boussinesq system 

In the case with periodic boundary conditions, a further simplifi- 
cation consists in restricting the model to those parts of the mean 
quantities which are independent of z. This is equivalent to re- 
defining the average as a volume rather than a horizontal one. 
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We obtain 



'R-xy — 

V — 



^yy 
7? 



xy 

-ouFxQu. 

4£l X TZy Z 

-OuFydz + 2fl x (TZ zz 

-2aT z g z - 4fl x TZ yz 



'-TZ 3 ' 2 , 



3L 

^QxT&xy H~ ^Qz'^yz ~ ■^TZ^'X 



A^l z ^'xy A-lZ^~yy 



3L ' 



, — 2fi 7? 



TZ 

a u tz 



(B.l) 



3L ' 



Ty = 

Tz = 
Q = 



-TZ xz Gq + 2^l z Ty — AjrT x , 



—lZ yz Go 



2^\ X J~ g < Z^1 Z J~ X AjzJ~ y, 

TLzzGq - aQg z - 2fl x T y - AjrJ" 2 , 
2T Z G - AqQ, 



(B.2) 



with Go = ATq/L z which is for ft = equivalent to Eqs. (53) 
of GOMS10. The resulting equation for TZ reads 



n 



-2aT z g z 



(B.3) 



and is not explicitly influenced by rotation. 

In the non-rotating case the equations for TZ XX , TZ yy , TZ ZZ , 
T z and Q form a closed system which can be solved in sepa- 
ration from the remaining equations. Once the solution of the 
former is known the latter can be solved where one finds again 
two separate systems: {IZxz^x} and {TZy Z , Ty}. The latter two 
systems have the same shape, and when assuming that there 
is a stationary solution for TZ from the first system, we ar- 
rive at the eigenvalue problem for the growth rate of an ansatz 
r R-xz,T x ,TZ yz ,F y ~ exp(Ai) 



An + A ag z 

G Ajr + X 



= 



(B.4) 



with constant An, Ajr. The solutions are 

An + Ajr l{A n -A T f 



Ai. 



± 



ag z G Q (B.5) 



and given that ag z Go > for convection, unstable solutions 
cannot completely be ruled out for sufficiently large values of 
this product, but had most likely to be considered unphysical. 

Nontrivial closed form stationary solutions can be derived 
for the special settings $1 = and $ = (pole): In both cases 
we have 



7? — 7? 



hence 



^-xx — 

TZ = 

Hzz = 
Q = 



C 2 K^ 
vv ~ 3L A n 



C 2 



2ag z G L 2 

CiCq 
3Ci + C 2 
3(d + C 2 ) 

-——■—7?. 

ag z C 7 



Ci 
C 7 



3(Ci + C 2 ) 



3Ci + C 2 



TZ 



TZ, 



3(d 
T z = —. 



C 2 ) 
2ag z L ' 



(B.6) 
(B.7) 

(B.8) 
(B.9) 
(B.10) 



which coincides with the solution given in GOMS 10. In turn it is 
under these conditions possible to determine the {Ci} uniquely 
when 1Z XX — TZ yy , TZ ZZ , T z and Q are given from a DNS: 



(B.ll) 



3TZ xx Ci — (1Z ZZ — 1Z XX ) C*2 — 
3TZ ZZ C 1 + {3TZ ZZ - TZ)C 2 = -Qag z LF z /n 1/2 
C 6 = ~(G K ZZ + ag z Q)L/T z n^ 2 
C 7 = -2G LT Z /Q1Z 1/2 . 
Inserting (IB. 81 ) in dB.5t we obtain 
(A K - A^) 2 + Aag z G a - (A n + Ajr) 2 = 4(ag z G Q - A n Ajr) 

C\ + c 2 c 2 

CV 3C7 



= 2ag z G 1 



(B.12) 



the sign of which depends solely on the parameters {G\} and 
not on ag z Go. Requiring (IB. 12b to be negative provides an ad- 
ditional constraint. A corresponding generalized condition, en- 
suring overall stability, is referred to in Sec. 13.2. ll 

Another special situation is found at the equator ($ = 0, 
hence fl z = 0) where the system ( IB. Il l decomposes into a closed 
one for the quantities 7Z XX , 1Z yy , 1Z yz , 1Z ZZ , T y , T z and Q and 
another one for 1Z xy , 1Z XZ and T x which can be solved once TZ 
from the first system is known. The latter reads in the stationary 
case 



= 2fl x 1Z xz — AnTZ xy , 
= -aF x g z - 2^l x 1Z xy 
= ~1Z xz Go — AjzT x , 



AnK 



n'^-xz, 



with T x being eliminated by the last line by J- x = 
— Go TZ X z j An- The remaining two equations form a homoge- 
neous linear system for !Z xy and 1Z XZ having the determinant 



ag z G Q 



Ch + Ch 



c x + c 2 



Nontrivial solutions would be possible if TZ were to assume a 
special value depending on the parameters. However, this had 
the unphysical consequence, that then TZ xy , TZ XZ and T x would 
depend on an arbitrary quantity. So we have to conclude, that 
they either vanish or are time-dependent. In the latter case we 
had to require stability, so these quantities were bound to decay 
to zero or to perform stationary oscillations with an arbitrary 
amplitude. As the only physically meaningful option we assume 
that they vanish. 

The remaining system reads 



AnTZ x 



3L /C ' 



4n x TZ yz - AnTZ. 



-aT y g z ■ 
-2aF z g z 



n'^yy 1- 

AQxT^yz 



9ln 5 ' 2 

3L ' 

— TZyy) — AnTZ yZ , 



A n n z 



3L 



TZ 3 / 2 , 



-TZy Z Go + 2il x J 7 z — AjzTy, 

HzzGq - aQg z - 2Vt x F y - AjzT z 
-2F z Go~A Q Q. 



From the first line it follows TZ = Ci 2 (TZ yy + TZ ZZ ), G\ 2 = 
3(Ci + C 2 )/(3Ci + 2C* 2 ) and from the last Q = -2G Q F Z /A Q 
leaving a system with five variables only. It can be broken down 
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to a nonlinear equation for TZ which is solved by TZ = and the 
solutions of 

2a x (2K-TZE{TZ)/C 12 )-A n KD(TZ)-ag z G(TZ)E(TZ) = 0, 

(B.13) 

completed by 



T z = KTZ 3/2 , T y = G(K)K 



3/2 



KD{n) 



_ TZ K „ /2 ^ 

VV ~C 12 E(TZ) ' yz E(1Z) 

TZ - K TZ 3/2 K-^-- Cl+ ° 2 

zz ~ E(1Z) ' 3L C 12 L ' 



with 

E(K) = A Z (K)D(K) + B Z (TZ), G{1Z) = + 



K AyD + By 



D(TZ) = ag z 
Go , 



4n x B z (1Z) + A n By(1Z) 



Wni + A 2 n - ag z (4n x A z (K) - A n A y {K)) 



A v = -f{^9zG /A Q - Aj?), B y = -G 2fl x /A = -A z , 



B, 



-GqA^/A, A = 4Ql+A^(A^-ag z 2G /A Q ) 



It cannot be guaranteed that (IB. 13b has positive solutions for TZ 
for any arbitrary set of parameters, in particular for arbitrary pos- 
itive {G t }. 

In contrast, for fl ^ and i3 ^ 0, 90° none of the compo- 
nents of TZ and T disappear and the determination of the {Cj} 
from DNS results has to deal with an overdetermined system: 10 
equations vs. 4 unknowns. 

With respect to the realizability constraint Q21 l. an analysis 
analogous to that of GOMS 10, Appendix A, but with rotation in- 
cluded, leads to the following relation for the temporal derivative 
of the quantity T = XfT^Xj = X^lZij - Q~ 1 J r iJ r j)X j 



DT = [T ■ X) 2 ^ L vlZ (2C 6 -Ch-d- C 2 ) 

+ ^ n 3/2 X 2_ {Ci + C2) TVn_ TdkJJk (B14) 

+2{T ■ X)Q- 1 X J T jk d k (e + T ) - IXiSukQtTjkXj 

where D = dt + V A. Repeating the arguments of GOMS 10 
here, one finds that the realizability condition is not affected by 
the presence of rotation, since in (IB. 141 ) is multiplied by the 
vanishing term TijXj. Similarly, by retaining the model coeffi- 
cients C v , C K and C„ K one can derive the following expression 



DT = [T-X 



2Cq — C7 — Ci — C 2 



2{v + n)C UK - kC k - vC v 



-^R?' 2 X 2 - (Ci + C 2 )I^ - Td k U k 

-2{T ■ X)Q- 1 X 3 T 3k d k (Q + To) - 4X: i e ilk il l T jk X j 

vG^ r 
L 2 



X{ T 9 T-ij Xj , 



from which one obtains the realizability criterion 



2Cq— Cj— C\— C 2 + 



2{u + k)C v 



KC K - VC V 



> 0. (B.15) 



Apart from being explicitly position-dependent via L(z), this 
criterion cannot be formulated as a condition for the model pa- 
rameters alone, unlike (132V However, we can infer the two suffi- 
cient conditions d32l i and 2(y + k)C„ k — nC K — vG v > 0. With 
Pr = 1 the latter one can be written as 4C„ K — G K — C v > 0, 
which is satisfied by the values G UK « 6, C v « 12 and C K w 2 
given in GOMS 10. 
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Table B. 1 : Summary of the Boussinesq DNS results. Normalizations (indicated by a tilde) are carried out with agdATo for Reynolds 
stresses, with ATo(agdATo) 1 / 2 for heat fluxes, and with (AT ) 2 for temperature variance. 
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83 


0.128 


0.738 


1.376 


0.131 


-2.688 


0.282 


-0.721 


-6.584 


20.468 


0.286 


D12 


45° 


10.24 


0.89 


91 


0.152 


0.816 


2.599 


0.152 


-2.936 


0.348 


-0.114 


-7.675 


24.670 


0.343 


D13 


45° 


12.96 


1.39 


66 


0.076 


0.677 


-1.131 


0.098 


-3.511 


0.169 


-1.734 


-5.371 


14.012 


0.190 


El 


60° 


0.04 


0.06 


87 


0.112 


0.042 


-0.168 


0.121 


-3.609 


0.363 


-0.190 


-3.910 


29.059 


0.344 


E2 


60° 


0.16 


0.13 


77 


0.091 


0.175 


-0.350 


0.114 


-4.614 


0.259 


-0.460 


-5.583 


21.095 


0.263 


E3 


60° 


0.36 


0.22 


69 


0.078 


0.399 


-0.549 


0.105 


-4.219 


0.194 


-0.791 


-5.760 


16.072 


0.211 


E4 


60° 


0.64 


0.32 


63 


0.068 


0.595 


-0.712 


0.097 


-3.530 


0.150 


-1.172 


-5.590 


12.584 


0.174 


E5 


60° 


1.00 


0.43 


59 


0.063 


0.690 


-0.781 


0.087 


-2.856 


0.128 


-1.355 


-5.060 


10.734 


0.154 


E6 


60° 


1.44 


0.53 


57 


0.060 


0.738 


-0.729 


0.081 


-2.384 


0.115 


-1.364 


-4.700 


9.379 


0.139 


E7 


60° 


1.96 


0.60 


59 


0.064 


0.748 


-0.565 


0.084 


-2.232 


0.123 


-1.330 


-4.817 


9.763 


0.146 


E8 


60° 


2.56 


0.67 


61 


0.070 


0.819 


-0.510 


0.087 


-2.015 


0.135 


-1.257 


-4.653 


10.379 


0.156 


E9 


60° 


4.00 


0.76 


67 


0.084 


0.781 


-0.243 


0.101 


-1.838 


0.170 


-0.851 


-4.983 


12.298 


0.185 


E10 


60° 


5.76 


0.82 


74 


0.101 


0.906 


-0.135 


0.118 


-1.880 


0.213 


-0.872 


-5.667 


15.029 


0.226 


Ell 


60° 


7.84 


0.87 


82 


0.120 


0.809 


0.087 


0.138 


-1.913 


0.268 


-0.317 


-6.481 


18.542 


0.278 


E12 


60° 


10.24 


0.90 


90 


0.140 


0.925 


-0.002 


0.160 


-2.037 


0.336 


-0.390 


-7.392 


22.716 


0.337 


E13 


60° 


12.96 


0.95 


96 


0.158 


0.784 


-0.150 


0.179 


-2.009 


0.397 


-0.083 


-8.190 


26.151 


0.390 
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Table B.2: Summary of the Boussinesq DNS results continued. For explanations see Tabl dB . 1 1 



Run 





Ta/10 6 


Co 


Re 


T^xx 


TZxy/W 2 


7^/10 2 


Ttyy 


Hyz/W 2 


Hzz 


•7V 10 


Ty/1Q 2 






Q 


Fl 


75 


0.04 


0.06 


85 


0.108 


0.015 


0.001 


0.120 


—3.963 


0.347 


—0.028 


—4.396 


27.848 


0. 


.331 


F2 


75 


0.16 


0.13 


75 


0.089 


0.089 


—0.175 


0.114 


—4.675 


0.246 


—0.187 


—5.806 


20.114 


0. 


.253 


F3 


r7rO 

75 


0.36 


0.23 


66 


0.071 


0.287 


—0.332 


0.104 


—4.106 


0.171 


—0.493 


—5.951 


14.343 


0. 


.193 


F4 


To 


0.64 


0.33 


60 


0.061 


0.366 


—0.314 


0.094 


—3.267 


0.134 


—0.575 


—5.607 


11.250 


0. 


160 


F5 


75 


1.00 


0.44 


57 


0.055 


0.467 


—0.366 


0.086 


—2.566 


0.115 


—0.674 


—5.068 


9.405 


0. 


140 


rO 


i 


1 AA 


u.oo 




U.UOl 




— U.oZ ( 


fl 078 
U.U ( o 


— i .yoo 


U. 1UU 


— u.ooy 


A A*\A 


7 Q1 Si 


u. 




F7 


75° 


1.96 


0.64 


55 


0.053 


0.659 


-0.399 


0.081 


-1.683 


0.107 


-0.701 


-4.287 


8.072 


0. 


.127 


F8 


75° 


2.56 


0.72 


56 


0.057 


0.812 


-0.516 


0.084 


-1.432 


0.111 


-0.770 


-4.133 


8.055 


0. 


130 


F9 


75° 


4.00 


0.84 


60 


0.067 


1.007 


-0.519 


0.093 


-1.075 


0.127 


-0.686 


-3.948 


8.733 


0. 


.144 


F10 


75° 


5.76 


0.89 


69 


0.085 


1.218 


-0.495 


0.112 


-1.113 


0.175 


-0.518 


-4.467 


11.654 


0. 


.187 


Fll 


75° 


7.84 


0.93 


76 


0.102 


1.241 


-0.473 


0.132 


-1.203 


0.227 


-0.229 


-5.079 


14.846 


0. 


.233 


F12 


75° 


10.24 


0.98 


83 


0.118 


1.241 


-0.693 


0.147 


-1.240 


0.279 


-0.171 


-5.637 


17.673 


0. 


.276 


F13 


75° 


12.96 


1.02 


89 


0.134 


1.402 


-0.606 


0.164 


-1.138 


0.332 


0.196 


-5.919 


20.774 


0. 


319 


Gl 


90 


0.04 


0.06 


85 


0.108 


0.015 


—0.025 


0.120 


—3.925 


0.347 


—0.046 


—4.319 


27.933 


0. 


.332 


G2 


91) 


0.16 


0.14 


75 


0.087 


0.010 


—0.064 


0.114 


—4.725 


0.240 


—0.073 


—5.925 


19.723 


0. 


.250 


G3 


91) 


0.36 


0.23 


66 


0.069 


—0.012 


0.004 


0.107 


—4.144 


0.168 


—0.008 


—6.277 


14.067 


0. 


.192 


G4 


90 


0.64 


0.34 


60 


0.058 


—0.010 


—0.002 


0.094 


—3.097 


0.130 


—0.018 


—5.476 


10.860 


0. 


155 


G5 


90 


1.00 


0.46 


55 


0.050 


—0.029 


0.045 


0.084 


—2.366 


0.108 


0.053 


—4.896 


8.771 


0. 


132 


G6 


90 


1.44 


0.57 


53 


0.045 


0.004 


0.018 


0.082 


— 1.870 


0.098 


0.028 


—4.503 


7.543 


0. 


119 


G7 


90° 


1.96 


0.68 


52 


0.040 


0.036 


-0.020 


0.085 


-1.381 


0.090 


-0.067 


-4.045 


6.579 


0. 


109 


G8 


90° 


2.56 


0.85 


48 


0.026 


0.015 


0.013 


0.086 


-0.753 


0.068 


-0.003 


-2.888 


4.573 


0. 


083 


G9 


90° 


4.00 


0.93 


54 


0.032 


-0.032 


0.054 


0.103 


-0.725 


0.098 


0.075 


-2.992 


5.823 


0. 


103 


G10 


90° 


5.76 


0.98 


62 


0.039 


-0.017 


0.032 


0.128 


-0.811 


0.135 


0.029 


-3.692 


7.680 


0. 


133 


Gil 


90° 


7.84 


1.09 


65 


0.038 


-0.009 


0.025 


0.143 


-0.640 


0.153 


0.010 


-3.738 


8.229 


0. 


.146 


G12 


90° 


10.24 


1.16 


70 


0.042 


-0.001 


-0.021 


0.161 


-0.677 


0.184 


0.059 


-4.100 


9.925 


0. 


174 


G13 


90° 


12.96 


1.22 


75 


0.048 


0.008 


-0.070 


0.173 


-0.806 


0.220 


-0.045 


-4.411 


11.906 


0. 


207 



